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Abstract. We propose a method for simulating the behaviour of small clusters 
of particles that explicitly accounts for all mean-field and binary-correlation 
effects. Our approach leads to a set of variational equations that can be used to 
study both the dynamics and thermodynamics of these clusters. As an illustration 
of this method, we explore the BCS-BEC crossover in the simple model of four 
fermions, interacting with finite-range potentials, in a harmonic potential. We 
find, in the crossover regime, that the particles prefer to occupy two distinct pair 
states as opposed to the one assumed by BCS theory. 
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1. Introduction 

Feshbach resonances are now widely used to control the interactions between atoms 
in ultra-cold atomic gases. At temperatures low enough for superfluidity to occur, 
these resonances allow experiments to realise the crossover between a supcrfluid of 
Cooper pairs (BCS) and a Bose- Einstein condensate (BEC) of molecules (O'Hara 
et al. 2002, Chin et al. 2004, Kinast et al. 2005, Bourdel et al. 2003, Regal 
et al. 2005; Zwicrlcin et al. 2005). Despite breathtaking advances in experiments, 
the task of formulating a theory justified in connecting the two extremes remains. 
The difliculty stems from the increasing importance of correlations within the system 
and requires a theory which accounts for many-body effects at a level well beyond 
mean-field theory. 

There have been a variety of proposals for treating these correlation effects mostly 

based on the methods of diagrammatic many-body perturbation theory and relying 
upon the inclusion of various subsets of interaction diagrams. Unfortunately, due to 
the lack of a clear expansion parameter, it is difficult to justify these approaches. 
Moreover, at T = 0, simple BCS mean-field theory smoothly interpolates throughout 
the BCS-BEC crossover despite it's neglect of correlations. This behaviour was first 
noted by Leggett (Leggett 1980), but the utility of BCS theory is apparent when 
compared to the results of recent experiments (Regal et al. 2005) and Monte-Carlo 
calculations (Astrakharchik et al. 2005) of the crossover. BCS theory continues to 
work surprisingly well within regimes where a mean-field theory is hardly justified 
(we will return to this point later) . At finite-temperatures it is true that conventional 
BCS theory diverges in the crossover, but a simple Gaussian expansion about the 
mean-field theory fixes this problem (Nozieres 1995). The result is we have a theory 
that smoothly interpolates throughout the crossover region despite a the neglect of 
quantum correlations. This surprising robustness of BCS theory has motivated us 
to reconsider the mean-field approach and to look for a method to extend it in a 
physically transparent way. 

In order to understand the difficulties that correlations introduce let us recollect 
why mean-field theory is such an attractive approach to a many-body problem. The 
allure of mean-field theory is that it lets us squeeze all the complications of the 
many-body problem into a single interaction term leaving us with an effective single- 
particle problem. For a dilute gas, in order to include the effects of interactions 
among the particles, the only microscopic input needed comes from an Tinderstanding 
of the binary scattering problem. Unfortunately, as higher-order correlations form 
throughout the system, i.e., beyond binary, mean-field theory breaks down, and a 
knowledge of two-body interactions is no longer sufficient. Rather, for a complete 
theory, we must account for an increasing number of N-body correlations. A recent 
example that clearly illustrates this need is the four-body calculation that was 
necessary to find the correct molecule-molecule scattering within a dilute gas of 
fermions far on the BEC side of the resonance (Pctrov et al. 2004) . 

Since interactions increase the number of N-body correlations that govern the 
microscopic physics, a starting point for understanding the strongly-correlated regime 
is to better understand the physics of small clusters of atoms. In this paper we 
will study a cluster of four particles, interacting with finite-range potentials, at all 
interaction strengths across a resonance. This forms the simplest model one could 
imagine of the BCS-BEC crossover. On the BEC side we will have two interacting 
bound molecules, while on the BCS side we have Pauli blocking of the lowest energy 
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level resulting in an analogous "Cooper pairing." We will attack this problem with 
an approach that, although approximate in its solution of the four-body problem, can 
be generalized to larger clusters of particles without much difficulty. In section |21 we 
discuss the basis of our approach to studying clusters, with the aid of a generalized 
pair wavefunction, and contrast this wavefunction with the more restrictive BCS 
wavefunction. We further motivate this generalized pair wavefunction by showing 
how it describes both the BCS and BEC limits. We present a variational method for 
working with such a wavefunction in section and discuss our choice of trial pair- 
orbitals. Finally, in section ^ we discuss our numerical results for four interacting 
fermions, in a harmonic trap at T = 0, throughout the crossover regime and compare 
this to the BCS results. We find that the atoms tend to form pairs of varying size 
within the crossover as compared to the single width of the BCS pairs. We then 
discuss the relevance of this result. 

2. Binary correlations and pairing wavefunctions 

2.1. Hartree-Fock theory and binary- correlations 

In developing a controlled method for treating correlations let us first review Hartree- 
Fock (HF) mean-field theory. In the spatial representation of HF theory, the HF 
wavefunction is an anti-symmetric function of N single particle orbitals: 

= A{0i(l)(/.2(2)...0^(iV)}, (1) 

where the operator A enforces the anti-symmetry of the overall wavefunction. A 
wavefunction formed from single-particle orbitals, as in equation when inserted in 
the Schrodinger equation, can only account for single-particle effects, i.e., the mean- 
field potential felt by each particle. Correlations, on the other hand, are multiple- 
particle effects and may be defined as all many-body effects neglected by the HF 
approach. We generalise equation to account for the next order, or all binary 
correlations, with an anti-symmetric function of N/2 — M two-particle orbitals or 
pair-orbitals: 

*P = A{0i(I, l')</>2(2, 2')... (/.m(M, M')}. (2) 

This wavefunction looks very similar to the BCS wavefunction. In its first quantised 
form, the BCS wavefunction, for a two spin mixture of fermions a = {], J,}, is formed 
by pairing each fermion with a partner of opposite spin, designating the pair by an 
orbital 0, and writing the full many-body wavefunction for the N fermions as the 
anti-symmetric product of M = N/2 pair-orbitals, i.e., 

*scs = ^{0(1, l')'^(2, 2') . . . 0(A^, M')}. (3) 

A remarkable feature of this wavefunction is that, despite the fact that the number of 
states available to each individual particle is severely restricted by Fermi statistics, a 
single pairwise orbital can become macroscopically occupied. It is this feature of the 
many-particle wavefunction that is at the origin of the coherent behaviour associated 
with superfluidity. 

We might, however, wonder if all the pairs will occupy just the one state. Indeed, 
as the inter-particle interactions are increased, and more states within the Fermi 
surface become available, we might well expect to generate pairs in more than one 
orbital. This additional freedom is already incorporated in the pair wavefunction 
of equation 10). If the BCS wavefunction is thought of as a HF wavefunction built 
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upon the Cooper pairs, the pair wavefunction may be thought to generate a "multi- 
configurational" BCS theory. The first question that we must now ask is how does 
one generate such a wavefunction? 

For our discussion, we will treat the case of a two-spin component gas of fermions. 
We will take pair states that are formed from a product of an even spatial wavefunction 
and an odd singlet spin-orbital. This assumption greatly simplifies the problem of 
anti-symmetrization. In fact, we are then able to write the pair wavefunction in the 
following determinantal form: 



(/.P(i)(l,l') 0P(i)(2,l') ... (/.p(i)(M,l') 
0P(2)(1,2') <^P(2)(2,2') ... (/.p(2)(M,2') 



•(4) 



0P(M)(1,M') cf>pi^M){2,M') ... cl)p(M){M,M') 

The sum in equation ^ denotes a summation over permutations of the orbitals 
1,2,... M. For example, in the case of 4-particles, the many-body pair wavefunction 
may be written as: 

* = 0i(l, l')02(2, 2') + 0i(2, 2')02(1, 1') - 0i(l, 2')'^2(2, 1') - <^i(2, l')02(l, 2'). (5) 

The pair wavefunction of equation ^ was explored extensively in the early days 
of superconductivity to explain the phenomena in terms of "Schafroth condensation," 
which is an analogy to Bose condensation for pairs of fermions. This wavefunction 
is also an excellent candidate for studying the BCS-BEC crossover since it can be 
shown to correctly tend to both the weakly interacting limits of the crossover without 
restricting the particles to occupy only a single orbital (Blatt 1964) . Unfortunately, 
equation ^ is rather difficult to work with since the number of determinants needed 
to generate the wavefunction scales like Ml. This difficulty will ultimately limit the 
number of particles that we can study with such a trial wavefunction. However, since 
we are only concerned with small clusters of Fermions, equation |5J) seems an excellent 
choice for a trial wavefunction. 

As mentioned above, the pair wavefunction is similar in form to the Hartree-Fock 
(HF) wavefunction. The HF wavefunction, however, is an anti-symmetrized product 
of single-particle orbitals and neglects all correlations apart from those due to the 
Pauli exclusion principle. In contrast, the pair wavefunction is an anti-symmetrized 
product of two-particle orbitals and may be roughly thought of as a second order HF 
wavefunction which accounts for the additional effect of all two-particle correlations. 
We now show how this trial wavefunction can treat both the BCS and BEC limits. 



2.2. BEC limit 

To illustrate the utility of equation (j^J, we shall examine the four-particle case, 
borrowing much of the following discussion from Blatt (Blatt 1964). To begin, we 
should note that the norm of the four-particle wavefunction has three contributions 
i.e., 

j \^\'^d^xd^y = A{N + D-S). (6) 

The first contribution, N — ((/)i|(/)i)(02|02), arises from the product of the norms of 
the pair-orbitals. The remaining two terms come from exchange processes inherent in 
the wavefunction. When both a spin up and a spin down are exchanged between pairs, 
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the resulting contribution to the norm is given by the term D = | (0i |(/'2) P which is the 
modulus square of the overlap between the two pairs. If only a single up or down spin 
is exchanged, however, we get the following contribution which we explicitly write: 

J </>t (1, 1>^(2, 2') (0i(l, 2')02(2, 1') + ^')Ml, 2')) d^xd'^y.i?) 

To understand the importance of these two exchange terms, let us consider the 
following. Assume for the moment that the integral of equation Q gives zero. The 
resulting norm of the many-body wavefunction is identical to the result we would have 
attained had we started with the much simpler wavefunction 

^b{1, 1'; 2, 2') « l')^2(2, 2') + M2, 2')M^, 1'). (8) 

Note, this wavefunction is now symmetric which means that we are in a regime where 
the fermion pairs behave like bosons. This shows that the contribution from S can 
be considered a measure of how much the pairs deviate from being bosons. The only 
place we would expect to find 5 = is in the BEC limit. The moment we move 
away from the asymptotic limit of BEC, this term begins to grow which means that 
a bosonic approximation to the pairing should in general be poor. Note, this also 
implies that two pairs do not scatter like two bosons, even when working with the 
BCS approximation to the wavefunction. 



2.3. BCS limit 

The pairing wavefunction of equation ||2Jl, which can describe the Bose condensed 
limit, is also appropriate in the limit of weak, attractive interactions where BCS theory 
should hold. This is clearly seen in the single orbital limit (i.e., 0i = (/)2 = . . .^p), 
equation (PJ, when the pair-orbital is given by 

0(x,z/)-^gfee''=l--^l. (9) 
fc 

For the non-interacting Fermi gas, the Fourier component represents a Fermi sea 
filled to Fermi wavenumber kf: 

3' = [o \k\ > kf ■ (10) 

The form of the pair-orbitals that arise in the superfluid state are virtually unchanged 
from the non-interacting case. The Fourier components of the orbitals only change 
for those with wave numbers close to the Fermi surface. This subtle difference results 
in a long-range tail to the pair-orbitals that correlate two fermions of opposite spin. 
These pairs have such a large average size that a great many of them may overlap 
which is in stark contrast to a dilute BEC of molecules. 

Despite the very different nature of a BCS superfluid from a BEC, we see that 
once again a wavefunction of the form of equation ^ is appropriate. As a flnal 
illustration of this point, we again refer to the case of four particles. Let us assume 
that, within the weakly- attractive limit, each pair-orbital asymptotes to a product of 
two single particle orbitals 

<Pi{x,y) 'fii{x)ipi'{y),(l)2{x,y) ifi2{x)(p2'{y)- (H) 
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Here we use the notation ipg- to denote single particle orbitals. 
equation may be written as 



Ml) ^i(2) 
¥^2(2) (/P2(2) 



9^1' (1') ^i'(2') 

ip2'{2') ^2' (2') 



In this limit, 



(12) 



which is again the result for the non-interacting case. 



3. Variational approach to the dynamics/thermodynamics 
3.1. Variational dynamical equations 

In the previous sections we have motivated the use of a pair wavefunction to describe 
superfluidity in both the strong and weak coupling limits. We must now address how 
such a wavefunction may be used to solve the many-body problem. In theory, one 
could self-consistently solve for the pair-orbitals and use the resulting wavefunction 
to calculate any desired expectation value. We will not perform such a calculation 
here, instead, we have chosen a simpler variational method of obtaining the pair- 
orbitals. We follow this approach for two reasons. First, a variational calculation is 
much simpler to implement than a full, self-consistent solution. Second, and more 
importantly, the variational approach we now describe will allow us to generalise our 
method to finite-temperatures. 

We begin our discussion by writing down the time-dependent Schrodinger 
equation in variational form 

f*^ d 

S^J^ dt{Q{t)\ih-~ H\Q{t)). (13) 

If we identify the kernel above as a Lagrangian, we are able to make a rather elegant 
connection with classical mechanics (for an excellent review of this subject see reference 
(Feldmeier & Schnack 2000)). Let us divide the Lagrangian in two parts as follows: 

C{Q{t),Q{t)) = {Q{t)\zhj^\Q{t)) - {Q{t)\H\Q{t)) 

= Co{Q{t),Q{t))-H{Q{t)). (14) 

The wavefunction is then parametrized in terms of the set of variational parameters 
{Q{t),Q{t)} = {qi{t),q2{t) . . . qN{t);qi{t),q2{t) . . . qN{t)}. Minimising equation (jTSJl 
results in the Euler-Lagrange equations of motion for the variational parameters 

ddC dC 

dt dqi dqi 

It is often more convenient to have the Eulcr-Lagrange equations in a form which only 
depends upon derivatives with respect to the parameters q^ and qa : 

^A^^g<, = -|^, (16) 

where A^^i is a skew symmetric matrix defined as 

= _ ^^^0 (17) 

dqudqa dqadqu 

The solution of equations (|15|l or (|16|l completely specifies the time-evolution of 
the quantum mechanical many-body problem. Besides providing a set of dynamical 
equations for the quantum mechanical system, these equations may be used to study 
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the system at finite-temperatures. By working within a microcanonical ensemble, 
assuming that the system is sufficiently ergodic, we may replace thermal averages 
with time averages. The evolution of the state \Q{t)) may be reconstructed from the 
variational parameters and then used to evaluate any desired expectation values. In 
this paper, however, we limit ourselves to the zero temperature, equilibrium case so 
will not concern ourselves with the time evolution of the parameters. 



3.2. Trial orbitals 

For the special case we are concerned with, that of a superfluid system, we will use the 
trial state \Q) = j^*), where the spatial representation of the ket I^E*) is given by Q. 
The pair-orbitals which construct our many-body wavefunction are fully parametrized 
by the variables q^- 

10) = 10(91,92... 9Ar)). (18) 

The efficacy of this approach is dependent upon a good choice of trial pair-orbitals. 
For our present study, we make a very simple choice for the pair-orbitals which is the 
set of unnormalised Gaussians of the form: 

(x,yfcW) = exp(tl|£-(i±|>!). (19) 

The variables aa{t) {Aa{t)) parameterise the relative (cm.) width of the pairs at 
positions x and y. We chose Gaussian functions so the wavefunctions become exact in 
both the non-interacting Bose and Fermi limits. Of course, a Gaussian decays much 
faster than the ipb ~ exp(— a;/T)/r form expected for a bound-state near threshold 
(Taylor 1972). The quantitative errors resulting from this choice of wavefunction will 
be greatest when we approach the resonance where the exponentially decaying tail of 
the Gaussian is unable to account for the actual spatial decay of the wavefunction. 
However, in the the strongly-interacting region, i.e., near resonance, we believe that 
Gaussian orbitals are again adequate. Although the binary scattering near threshold 
predicts a long tail on the wavefunction, the only length scale that should effect the 
many-body system is the inter-particle spacing. Since our results connect smoothly 
between the unitary and weakly-interacting regimes we believe that our choice for the 
orbitals will be able to display qualitatively accurate features of the crossover. 



4. Four-particle BCS-BEC crossover 
4.1. Four-particle model 

We performed a numerical minimisation of the variational energy 

{mm. (on) 
^^7*^' ^^^^ 

where the four-particle Hamiltonian is defined as 

^ = E (-£(^x. + K) + + y')) + E - y^D- (21) 

The trapping frequency is given by lj and we only allow for interactions between spin 

t and I particles. We made use of a Gaussian potential of the form: V{\Ki — yj\) — 

I I ^ / 

_[/e-|x.-y,i ^ where U gives the depth of the potential and a represents the width. 
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A Gaussian is chosen simply because it separates for each dimension, a fact that will 
greatly aid in the numerics. To tune the interaction we set the width cr to a given 
value and adjust the depth U of the potential about it's first bound state which 
generates an effective threshold scattering resonance. The s-wave scattering length 
has the dispersive form Gsc — — 7r^/^{/(T'^/(4(l — IJa^/2)), and we have absorbed the 
remaining physical constants into U = {rn/h^)U . 

The minimisation is straightforward to perform using a multi-dimensional quasi- 
Newton method (Press et al. 2002). The greatest difficulty posed by equation if^ . 
however, is in the evaluation of the matrix elements. Each integral must be evaluated 
over four particles as opposed to two which effectively doubles the dimensionality. 
Various exchange and spatial symmetries help to reduce the dimensionality of the 
problem, nonetheless, the matrix elements remain the most inefhcient part of our 
program. This has imposed two constraints: 1.) We have to work with a rather large 
width to the potential so we can integrate over enough points within the range of 
the potential. This, of course, degrades the universality of the potentials we have 
chosen. II.) The width of the orbitals cannot get much smaller than the width of the 
potentials or again we will not be able to accurately evaluate the integrals. Ultimately, 
this means that the distance we can move below the resonance is limited by our grid 
spacing. Both of these problems can perhaps be circumvented by a more efhcient 
evaluation of the matrix elements, e.g., through the use of a non- linear grid or an 
adaptive form of Gaussian-quadrature. Our current approach should, however, be 
sufficient to justify our qualitative results. 

We first confirmed that the orbitals of equation (|I9|I result in SItlo for the non- 
interacting ground state energy. By dropping terms related to the exchange of a single 
fermion between the pair of orbitals, we could also reach the non-interacting Bose limit 
Qhw (as explained in section [^J. After performing these baseline tests, we fixed the 
width of the interaction potential at a = 0.3ao, where ao = \/h/muj is the oscillator 
length. We chose this value of a as the minimum value for which we could accurately 
perform the matrix elements. 

^.2. Numerical results 

In the following analysis, we refer to the results obtained with the pair wavefunction 
of equation ^ as the "pairing model" while the "BCS model" refers to the use 
of equation JSJ. For the four-particle problem, the crucial difference between the 
two models is that, while the BCS model restricts the pairs to a single orbital, the 
pairing model allows for the occupation of two separate orbitals. Before discussing our 
numerical results, we should make two points. First, the variational orbitals we have 
chosen will modify the zero of the resonance, i.e., the point at which a two-body bound 
state first appears. To account for this shift wc calculate the potential depth where the 
first bound state arises, and then simply shift the zero of the resonance accordingly. 
The scattering length is then given by the theoretical value a^c about this point. For 
the Gaussian orbitals we have chosen, of range a = O.Saq, the bound state appears at 
a depth of Ub ^ 3.51 ftw. Note, we work in terms of the dimensionless length a/a^ as 
opposed to kfa. For four particles the Fermi vector kf is somewhat arbitrary whereas 
oo seems the more natural unit of length. Second, our results clearly show that, as we 
move toward weak-repulsive interactions, the four-particle cluster collapses. This will 
prevent us from reaching the dilute BEC regime, a result which was perhaps expected 
given the large spatial extent of our inter-particle potentials. One could argue that this 
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Figure 1. The main graph plots the total energy minus the binding energy of 
each pair e;, as a function of the inverse scattering length. The BCS result {□) 
is plotted with that of the two-orbital, pairing model (o). The insert shows the 
relative difference between the two energies. All energies are in units of the trap 
energy ttw and lengths are in units of the oscillator length oq . 



collapse behaviour is quite physical and that the experimentally observed crossover 
to a dilute BEC is only quasi-stable, the four particles eventually solidify and falling 
from the trap. We could perhaps simulate this quasi-stability by using very short 
range potentials so that the minimisation algorithm has a difficult time finding the 
true ground state but rather leads us to the local minimum of a quasi-stable, dilute 
BEC of molecules. 

In figure we compare the ground state energy of the pairing model, after 
subtracting off the variational bound state energy of the pairs, to the same quantity 
as given by the BCS model . As the resonance is approached from above, meaning we 
begin on the BCS side of the crossover, both models initially predict the same ground 
state energy. The pairing model, however, branches off toward lower energies, only to 
converge again with the BCS model well below the resonance. The insert of figure 
plots the difference between these two curves to emphasise their deviation which is 
notably peaked just below the resonance. It should be noted that both curves display 
a similar, linear behaviour in the universal regime near the resonance. Unfortunately, 
the collapse of the cluster, which drags all our curves downward as we move toward the 
BEC side of the resonance, prevents us from making a qualitative comparison between 
our model and other predictions for the equation of state within the unitarity limit. 
The energy continues to shoot toward negative values below the resonance signalling 
a complete collapse of the cluster. 

To better visualise the four-particle wavefunction, in figure (|5J) we plot the relative 
and cm. widths as a function of inverse scattering length. The relative width serves 
as a measure of the size of each pair. If we move toward the resonance from the BCS 
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Figure 2. The graph on the left shows the relative width of the pairs, given by 
the variational parameter o, as a function of inverse scattering length. The BCS 
result (*) is fitted with a line to guide the eye. The widths of the two orbitals of 
the pairing wavefunction are given by the symbols □ and o. The graph on the right 
shows the cm. width of the pairs which is given by the variational parameter A. 
The symbols arc the same as in the leftmost graph (note, the circles and squares 
overlap). Again, all lengths are in units of the oscillator length ao. 



side, one orbital is preferred by both pairs of atoms, which is rather large and stretches 
across the trap. As we continue to approach the resonance, the pairs separate into two 
orbitals with one larger and one smaller than the BCS orbital. Both orbitals gradually 
shrink as we move through the resonance only to converge well below the resonance 
at which point they are on the order of the size of the potential. We also see from the 
figures that the cm. widths of the pairs, in the pairing model, tend to spread more 
within the crossover region than does the BCS orbital which is what we might expect 
to accompany the diminishing Fermi sea. Both pairs, however, tend to lock onto the 
same cm. width throughout the crossover. The cm. width quickly tends toward zero 
as we move below the resonance which is another sign that our cluster is collapsing. 

5. Conclusions cind Outlook 

We have developed a variational method for studying small clusters of correlated 
fermions that accounts for all binary-correlations among the particles. The method 
is based upon a very general form of the many-particle wavefunction, similar to BCS, 
but not as restrictive in that it allows each pair to occupy a separate orbital. We have 
described in detail a variational method that can be used to study the equilibrium, 
dynamical and finite-temperature behaviour of these clusters using this wavefunction 
as an ansatz. As an example, we have applied our method to the study of the 
equilibrium properties of the BCS-BEC crossover for four particles, interacting through 
finite-range potentials, in a harmonic trap. We see signatures of universal behaviour 
for our four-particle system near the resonance, but are restricted from moving into 
the dilute BEC regime by the collapse of the cluster. We have found that the system 
fragments from the BCS prediction of a single pair to two pairs, one larger and the 
other smaller than the BCS pair. The presence of distinct pairs at T = should 
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be contrasted with the idea of "preformed" or "noncondensed" pairs that is often 
discussed in relation to the crossover. These pairs are often thought to arise from 
finite-temperature effects or fluctuations on the order parameter. Here, however, the 
distinct pairs that form within the crossover region are truly a ground state property 
of the system at zero temperature. 

There are a few final points we should emphasise. The behaviour of the pairs, as 
displayed in figure (0), hints at why BCS theory interpolates throughout the BCS-BEC 
crossover so well since the BCS result tends to act as an average of the various pairs 
that form. However, one would need to look at larger clusters to see if this pattern 
holds up. It is clear, nonetheless, that by allowing the pairs to occupy different orbitals 
one can better model the true ground state. Although the full pairing wavefunction of 
equation is much too complicated to work with for large numbers of particles, the 
inclusion of a restricted set of orbitals could be implemented as the variational state 
within quantum Monte-Carlo (MC) calculations. It would be of interest to see, for 
instance, the effect of an additional pair-orbital on the fixed-node MC results when 
a BCS wavefunction is chosen as the variational state determining the location of 
the nodes. As a final note, we feel that the results of the present work, due to the 
rather large width of the potentials, might be ideally suited for studying systems such 
as liquid He^ where the range of the potential is on the order of the inter-particle 
spacing. In these systems, the collapse we see of the clusters may be unavoidable and 
therefore of more interest than the quasi-stable, crossover state that is the current 
concern in dilute gases. 
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